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Global Superfluid Phase Diagram of Three Component 
Fermions with Magnetic Ordering 
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We investigate a three component fermion mixture in the presence of weak attractive interactions. 
We use a combination of the equation of motion and the Gaussian variational mean-field approaches, 
which both allow for simultaneous superfluid and magnetic ordering in an unbiased way, and capture 
the interplay between the two order parameters. This interplay significantly modifies the phase 
diagram, especially the superfluid-normal phase boundaries. In the close vicinity of the critical 
temperature and for small chemical potential imbalances, strong particle-hole symmetry breaking 
leads to a phase diagram similar to the one predicted by Cherng et al. [Phys. Rev. Lett. 99, 
130406 (2007)], however, the overall phase diagram is markedly different: new chemical potential- 
driven first and second order transitions and triple points emerge as well as more exotic second order 
multicritical points, and bicritical lines with 0(2, 2) symmetry. We identify the terms which are 
necessary to capture this complex phase diagram in a Ginzburg-Landau approach, and determine 
the corresponding coefficients. 

PACS numbers: 37.10.De, 74.25. Dw, 67.60.-g 



I. INTRODUCTION 

Experiments with ultracold atoms opened a fascinat- 
ing way to study strong correlations and the emer- 
gence of exotic phases in a controlled wayi Paradig- 
matic solid state physics models such as the fermionic and 
bosonic Hubbard models have been realized, Mott insu- 
lating and magnetic phases^ 3 - as well as various kinds 
of fermionio^— and bosonio 2 -^ superfluid phases have 
been observed. Topological excitations, e.g. vortices^, 
solitonsii, 2D and 3D skyrmionic excitations^, and knot 
configurations^ 3 - have been subjects to intensive research. 
Introduction of artificial gauge fields has also been con- 
sidered both theoretically and experimentally", indicat- 
ing that the realization of the quantum-Hall effect and 
related phenomena with cold atoms are within reach. 

Cold atomic systems provide, however, not only a 
way to study models emerging in solid state physics, 
but they were also proposed to be used to mimic phe- 
nomena appearing in high energy and particle physics. 
In particular, attractive three component mixtures have 
been proposed to simulate quark color superfluidity 15 
and "baryon" formation,— two fundamental concepts of 
quantum chromodynamics (QCD). An experimental re- 
alization of these mixtures is very difficult, but not hope- 
less: although three component systems are plagued by 
3-particle losses^ - — nevertheless, Fermi degeneracy has 
been reached in 6 Li systems^ which may be just sta- 
ble enough to reach interesting phases such as the trionic 
("baryonic") regime Ji Also, systems with closed s-shells, 
similar to Yb may provide an alternative and more sta- 
ble way to realize almost perfectly SU(N) symmetrical 
states^ii^ 2 

In this paper, we focus on the weak coupling regime 
of an attractive three component mixture, and study 
its low temperature color superfluid phases. Our main 



purpose is to study the effect of chemical potential dif- 
ferences, and provide a complete phase diagram for the 
SU(3) symmetrical interaction, which can be considered 
as the three component analogue of the famous phase 
diagram of Sarma^ Surprisingly, although several stud- 
ies have been reported so far, such a phase diagram has 
not been discussed in sufficient detail so far, not even 
in the weak coupling regime considered here. The first 
analysis of Ref. |l5| assumed complete SU(N) symmetry 
and has not considered the effect of different chemical 
potentials. It neglected furthermore the coupling be- 
tween ferromagnetic and superconduct ing order parame- 
tcrs. However, as later noticed in Refs.R6landl24l. SU(S) 
symmetry allows for a coupling between magnetic and 
superfluid order parameters, and the onset of superflu- 
idity is therefore naturally accompanied by a ferromag- 
netic polarizatio n 24 i 25 and possibly domain formation.— 
The consequences of such coupling have been explored in 
Ref. [2J in the immediate vicinity of the SU(3) symmet- 
ric phase transition using a Ginzburg-Landau approach, 
however, the regime of lower temperatures has not been 
investigated. 

Throughout this paper, we shall proceed in the spirit of 
local density approximation and focus on a homogeneous 
system of three interacting fermion species, described by 
the Hamiltonian 

H = J2 /"d 3 r*t Q ( r )( Wo _ Ma )^ a ( r ) (i) 

- J2 X »f3 /d 3 r^(r)*t (r )^( r )tf a(r ) _ 

Here \&£,(r) creates a fermion in a hyperfme state a = 
1,2,3 with corresponding chemical potentials, /i Q . The 
interaction between the species is assumed to be local 
and attractive (A a ^ > 0)^ Furthermore, throughout 



most of this work, we shall also assume SU(3) symmet- 
rical interactions, \ a ^p = A. This assumption is a valid 
approximation for the 6 Li system in the high magnetic 
field limit j22 and it would be certainly justified for Yb-like 
closed s-shell systems (but with attractive interactions). 
This assumption is certainly justified for Yb-like closed 
s-shcll systems, and is also a valid approximation for the 
6 Li system in the high magnetic field limits Although 
the scattering lengths in the lowest three hyperfine states 
are slightly different in the latter system, one can use 
radio frequency and microwave fields to make them equal 
up to ~ 0.1% accuracy^ 

The particular form of the single particle operator "Ho 
in Eq. ijl]) is not very important, since Ho enters the 
mean-field calculations only through the corresponding 
single particle density of states (DOS), for which we as- 
sume a simple form, p(£) = po (1 + 7 £) and a rigid band- 
width cut-off at £ = ±W. Keeping the linear term po7 £ 
is crucial: this term is the primary source of the coupling 
between ferromagnetic and superfluid order parameters. 
Note that in the small coupling regime only the DOS pp 
at the Fermi energy and its first derivative are expected 
to have considerable impact on the phase diagram, and 
therefore we do not need to go beyond this simple linear 
approximation. We should remark though that the inter- 
actions rcnormalize the chemical potentials, and there- 
fore the position of the renormalized Fermi energy, £p 
and the corresponding single particle density of states, 
Pp must be determined self consistently^ 

Although we also discuss to a certain extent the role 
of fluctuations in Section |Vj the bulk of this work con- 
sists of a mean-field analysis. Even this is, however, 
not entirely trivial. In the Hubbard-Stratonovich ap- 
proach of Rcfs. [3(J and [24| the decoupling of the interac- 
tion into ferromagnetic and superfluid parts suffers from 
a certain degree of arbitrariness.— Treating the ferro- 
magnetic and superfluid order parameters at equal foot- 
ing therefore requires care. Furthermore, at lower tem- 
peratures the second order transitions turn into first or- 
der transitions, and the free energy develops several in- 
equivalent local minima. To cope with these difficul- 
ties, we applied two complementary methods: an equa- 
tion of motion method, where vertex corrections are sys- 
tematically neglected, and a Gaussian variational ap- 
proach. Both approaches are exempt from the arbi- 
trariness of the Hubbard-Stratonovich transformation, 
account for the interplay between ferromagnetic and su- 
perfluid order, and, remarkably, they both result in the 
same self-consistency equations. However, the Gaussian 
variational approach goes beyond the equation of mo- 
tion method in that it also provides an estimate for the 
mean field free energy, and allows us to locate first order 
transitions. Since previous works indicate that the Fuldc- 
Ferrell-Larkin-Ovchinnikov (FFLO) phase with spatially 
varying order paramete r 32 ' 33 appears only in a tiny re- 
gion of the phase diagram,^ here we restrict our investi- 
gation to spatially homogeneous phases. We shall neither 
consider Breached Pair (BP) or Sarma phases ; 23 ' 34 since 
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FIG. 1. (Color online) Structure of the phase diagram for 
a constant DOS. Top: SF forms between fermions with the 
closest chemical potentials, whereas at higher differences the 
normal state (N) is favored. Bottom: Cut of the SF phase 
diagram for T < T c Apart from the special points (full and 
empty circles) SF order always forms in one of the channels 
(12), (23) or (31). At the point p x = p y = (full circle) 
the Hamiltonian is 5*17(3) symmetric, and the transition is 
described by an 0(6) critical point. The first order lines, sep- 
arating different SF phases terminate in second order critical 
points with 0(2, 2) symmetry (empty circles). 



these would require fermions of very different masses. 

Before we turn to the more detailed presentation of 
the calculations, let us summarize here our most im- 
portant results. In the small coupling limit, T c <C W, 
the phase diagram is expected to become universal for 
SU(3) symmetrical interactions: it should depend only 
on the dimensional temperature, T/T Cl the dimension- 
less chemical potential shifts, 6p a /T c , and the dimension- 
less particle-hole symmetry breaking, 7 = 7T C , defined in 
terms of the critical temperature T c at the SU(3) sym- 
metrical point, p a = p. Fig. [T] shows the corresponding 
schematic phase diagram in case of a particle-hole sym- 
metrical situation, 7 = 0. The bottom figure shows a 
finite temperature cut of the phase diagram as a func- 
tion of the chemical potential differences, 



Mx = (Mi - M2) 

Py = (P! + p 2 - 2p 3 )/VG , 

for a temperature T fixed somewhat below the SU(3) 




FIG. 2. (Color online) Numerically computed SF-N phase 
boundary as a function of the chemical potentials at con- 
stant DOS for \pF = 0.1, and 7 = £f = 0, corresponding 
to Tc U(:i) /W = 0.0076. The SF-N transition becomes of first 



order below a temperature T ' 



(horizontal dashed lines). 



Vertical lines denote the 0(2, 2) critical points of second (solid 
line) and of first order (dashed line). 
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symmetrical transition temperature, T c . In the various 
gray regions two species of the smallest chemical poten- 
tial difference pair up to form a superfluid (SF) state, 
while the third species remain gapless. This explains 
the star-like structure of the phase diagram: superfluid 
phases appear around regions, where two of the chem- 
ical potentials become equal. As we discuss later, the 
high ("hexagonal") symmetry of the figure is a direct 
fingerprint of the SU(3) symmetrical interaction, and a 
discrete particle hole symmetry. The superfluid state 
is destroyed, once all chemical potential differences be- 
come large compared to the condensation energy (white 
region). Close to T c the chemical potential driven SF- 
normal transitions are of second order (black lines) , just 
as in case of a two component mixture^ The transition 
between different SF phases is, however, always of first 
order (dashed lines). 

The phase diagram also exhibits some interesting 
points of special symmetry At the point \x x = fi y = 
the Hamiltonian is SU(3) symmetrical, and correspond- 
ingly, the phase transition at T = T c and \x x = jj, y = is 
described by an 0(6) theory (the six components corre- 
sponding to the real and imaginary parts of the superfluid 
order parameters). In three dimensions, this symmetry is 
spontaneously broken for T <T C . On the other hand, at 
the points indicated by white circles in Fig. [TJ the com- 
petition of two order parameters most likely leads to a 
so-called 0(2,2) critical behavior (see Section [Vj . 

Fig.[2]shows the numerically obtained phase diagram in 



FIG. 3. (Color online) Phase diagram in the vicinity of the 
SU(3) symmetric point in the absence of particle-hole sym- 
metry, 7/O. The sixfold symmetry of the phase diagram 
is destroyed. A higher DOS can make SF ordering favor- 
able in a channel not of the smallest chemical potential dif- 
ference, due to the gain in condensation energy. For the 
absolute values of the SF order parameters Ay see color 
code (right). Parameters at the SU(3) symmetric point: 
\ a ^ BPF = 0.112, jW = 0.5, T c /W = 0.011, fr/W = 0.24 
(half-filling). 



a 3-dimcnsional plot under the assumption of SU(3) sym- 
metric interaction and particle- hole symmetry (7 = 0). 
The dome- like structures correspond to superfluid phases 
with pairing in the (12), (23), and (31) channels. Below 
the horizontal dashed lines the chemical potential driven 
phase transitions become of first order, while above these 
lines they are of second order. These lines are thus the 
analogues of the critical point identified by Sarma.— The 
SF-normal transitions on the "roofs" of the domes belong 
to the 0(2) universality class, while the black solid lines 
correspond to 0(2, 2) critical points. Finally, the cross- 
ing of the black lines at (j, x = \i y = corresponds to an 
0(6) critical point. 

This rich phase diagram is further complicated if one 
allows for particle-hole symmetry braking, 7 7^ 0i^£ On a 
larger scale, the 7 7^ phase diagram looks quite similar 
to the 7 = phase diagrams, presented in Figs. [T] and [21 



however, the structure of the phase diagram changes in 
the close vicinity of the SU(3) symmetrical point. This 
is demonstrated in Fig. G2 where the central region of the 
phase diagram is shown for T = T c and T = 0.5 T c . The 
absence of particle- hole symmetry destroys the hexagonal 
symmetry of the phase diagram, and leads to a trigonal 
structure, as predicted by Chcrng et al.^ In this central 
region a higher DOS, — and thus gain in condensation 
energy — may make SF ordering favorable in a channel 
not of the smallest chemical potential difference. This 
effect is most spectacular at T = T c , where by shifting the 
Fermi energy of two species one can increase the critical 
temperature, and induce superfluidity (see Fig. [3l top) . 
We remark, however, that in spite of the relatively large 
particle-hole asymmetry introduced, this central region 
is typically quite small compared to the rest of the phase 
diagram, at least for weak couplings, T c <C W. The 
orientation of the phases is, however, opposite to the one 
predicted in Ref. l24t to obtain the same orientation, we 
need to flip the sign of the slope of the DOS, and assume 
a hole-like Fermi surface, 7 < 0. We must also add here 
that the Ginzburg-Landau action of Ref. [2J is unable to 
capture the endpoints of the "trigonal" region, and one 
must retain higher order terms in the action to account 
for these (see Section HV)) . 

The rest of the paper is organized as follows: In Sec- 
tion UH we introduce our mean-field methods. We also 
discuss the symmetries of the order parameters, leading 
to rather strong constraints on the form of the phase di- 
agram. In Section inil we present our main results on the 
SF phase diagram, with and without particle-hole sym- 
metry, and compare our findings with results on two com- 
ponent systems. In Section HVl we present the numerical 
Ginzburg-Landau expansion of the free energy around 
the SU(3) symmetric point, and identify the terms re- 
sponsible for the main features of the central part of the 
phase diagram. In Section [V] we discuss the effect of 
fluctuations in the special 0(2,2) symmetric bicritical 
points. In Section IVII we comment on the experimental 
rcalizability of an SU(3) symmetric system. Some of the 
technical details of our calculations can be found in the 
Appendices. 



II. MEAN-FIELD CALCULATIONS 



In this section, wc first use an imaginary time equation 
of motion (EOM) method to derive the self-consistency 
equations for the SF and magnetic order parameters. 
Then, to address the low temperature regime, where 
these equations have multiple solutions,— we also develop 
a Gaussian variational approximation. This approach 
provides an estimate for the free energy and enables one 
to locate first order transitions. 
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FIG. 4. (Color online) Omission of the vertex corrections 
in the connected four point functions (l.h.s.). Heavy lines 
denote the full propagators, and the square stands for the 
vertex contribution. 



A. Equation of motion technique 

To simplify our notation, let us first introduce the 6 
component Nambu spinor field 

$(x) = (*(x),&(x)) T ■ 

Here we used the compact notation x = (r, r) for the 
space and imaginary time coordinates. The correspond- 
ing 6x6 propagator matrix D(xi,x 2 ) = — (T r $(xi) o 
$t(x 2 )) contains the normal as well as the anomalous 
Green's functions of the fields vE , Q (x) and, assuming spa- 
tial homogeneity, also obeys D(xi,X2) = D(xi — x 2 ). 
In order to derive equation of motion for the propaga- 
tors, we start from the imaginary time equation of motion 
(EOM) of the fields, 



d T * a (x) = [H,* a (x)} 



(2) 



The EOM of the part -(T T ^ a ,(xi)vErt(x 2 )) of the prop- 
agator follows from Eq. ([2]), and reads 

{d Tl + Ho(vi) - fi a ) (T T * a (xi)*J(x 2 )) = S aP 6 XlX2 (3) 
+ J2 2A aT (T T #t (xi)# 7 (xi)# a (xi)*t (x 2 )>, 



with S XlX2 denoting the four dimensional Dirac-dclta 
function. Similar equations hold for the anomalous prop- 
agators, -(T T * a (xi)*^(x 2 )) and -(T T *t ! (xi)^(x 2 )). 
To make further progress, we simplify the four point 
functions appearing in these EOMs, by simply neglect- 
ing the vertex contribution, as shown in Fig. |4j This 
approximation is almost equivalent to the usual BCS ap- 
proximation, however, it goes beyond that, since it allows 
for the simultaneous appearance of different kinds of or- 
der parameters in an unbiased way. Furthermore, even 
in the simple SU(2) case, it also incorporates, e.g., the 
rcnormalization of the Pauli susceptibility at the mean- 
field level (see Section IIII Cp . With this approximation, 
the equation of motion become solvable, and the Nambu 
propagator is found to take the following form in Fourier 
space 



D(iw n ,k)~ 



iuj n - B(£ k ) 



(4) 



Here w n = (2n + 1) 7r T are fermionic Matsubara frequen- 
cies and the matrix B(£) is defined as 



, , _ (Z - A 2A 
B( ^= I 2A+ -(£-A* 



(5) 



The matrices 



A Q/ g = \ a jid a 



f>, (6) 

A Q/ g = (fi a + 2 2_^ A Q7 n 77 )<5 Q , / 3 — 2A„ / gn* /3 (7) 

7 

denote the SF order parameter— and the renormalizcd 
chemical potential, respectively. They arc defined in 
terms of the mat 
lous densities d, 



The matrices n and d can be used to describe magnetic 
and SF ordering, respectively. However, it is more nat- 
ural to use A and A as order parameters. Note that, 
according to Eq. , magnetic ordering implies a shift in 
the renormalizcd chemical potential A a/ g, and this shift 
can thus also be considered as a magnetic order param- 
eter. 

The expectation values Eqs. © and © are given by 
the propagator D(xi — X2) at equal times and equal po- 
sitions, and are thus determined by Eq. Q. Taking the 
inverse of Eq. @ and performing the Matsubara sum- 
mation over the frequencies ui n we obtain 



d£p(0/(B(0) 



u' 



n 
-d H 



— n 



-d 



(10) 



where p(£) denotes the DOS of Ho, and / stands for the 
Fermi function. Eqs. (J3][7J) and (fTU|) thus determine self- 
consistent ly the order parameters A and A. We solve 
these equations iteratively, starting from random initial 
conditions, and performing the integrals in Eq. (|10l) nu- 
merically. Notice that /(B(£)) is a matrix function, 
therefore, its evaluation requires numerical diagonaliza- 
tion of the Hermitian matrix B(£) for each value of £. 

We remark that the matrix B(£) in Eq. ((SJ possesses 
a symplectic symmetry 



1 

1 



B(0 



1 

1 



= -B T (0 



(11) 



since the order parameters A and A are skew-symmetric 
and Hermitian, respectively. This symmetry makes the 
eigenvalues of B(£) come in pairs, (v, —v), and thus sim- 
plifies some of our calculations of the free energy in the 
next subsection. It is also responsible for the structure 
of the equal time, equal position propagator in Eq. (JTT 



B. Gaussian variational approach 

To investigate the low temperature phase diagram, we 
employ a variational method. This method consists of 
finding the best Gaussian approximation to the free en- 
ergy F = — TlogZ of the system. As a first step, we 



express the grand canonical partition function Z as a 
functional integral 



Dip Dip e 



-Sbl>,4>] 



(12) 



with the action written as S = So + Sj n t, and the non- 
interacting and interacting parts defined as 



nx 01 densities n, ana zn&z 01 z 


ne anoma- 


So 


n a p s^tfci)*^*!)), 


(8) 


Sjnt 


d a/3 =(*„(aJi)*/j(a;i)). 


(9) 





\l 



dld2^(l)P o - 1 (l,2)0(2), 



(13) 



^A Q/3 <ixip a (x)ipp(x)ip (x)ip a (x). (14) 

a/3 J 



Here <j> = (ip, ip) is a Nambu spinor field and we used 
the notations "1" = (ri,n,vi), and Jdl..., to denote 
the integration over space and imaginary time variables 
and the summation over Nambu indices (v\ = 1, ... ,6) 
in a compact way. The inverse propagator 



-^0 — Oxixt I CV2 



Ho- fl 






-(Ho -ft) 



(15) 



contains the single particle Hamiltonian of the free fields, 
Ho, where p, a p = fi a 8 a p is a 3 x 3 diagonal matrix con- 
taining the chemical potentials. 

Our Gaussian approximation of the free energy is based 
on the standard inequality 37 



F < F G [V] = -TlogZ v + T(S - S v )v 



(16) 



Here the partition function Zt> and the average (. . . }x> 
are defined in terms of the Gaussian action 



S V = -\J Ald2<p(\)V~ 1 (\,2)<p(2). (1 



Z-d = / Dip Dip e 



-s-DbP,i>] 



1 



)v ~Z~^' ^^ ' 



S-dWM 



7) 
(18) 
(19) 



Since we do not want to restrict our investigations to 
actions that can be associated with a Hamiltonian, we do 
not require Sp to be local. Nevertheless, at the saddle 
points of Fg, St> turns out to be local, and there exists a 
Hamiltonian associated with it (see Eqs. (|22j|24p below). 

Since the action Sx> is quadratic, the propagator ma- 
trix of the Nambu fields can be written as 



V(l,2) = -(<P(l)<f>(2)h 



(20) 



and expectation values can be evaluated using Wick's 
theorem. We remark that the choice (|2TJ]) automatically 
fixes a certain ambiguity in the definition of D^ 1 . (For 
details see Appendix [Cj) The best Gaussian approxima- 
tion is given by the minimum of the functional Fq\V], 
where Fg satisfies the saddle point equation 



SFr, 



5V(1,2) 



0. 



(21) 



As is shown in Appendix [C] this equation is equivalent 
to the self-consistency equations (J5I6I7P and (fTU)) of the 
EOM technique, and amounts in V~ l being a local, 



V- 1 (l,2) = 5(x 1 ~x 2 )T>- 1 (x 2 ). 



(22) 



with the matrix operator on the r.h.s being just the in- 
verse propagator Eq. (|4]) in real space, 



D 



Or 



Ho(r 2 )-A 2A 

2A+ -(U (r 2 )-A*) 



(23) 



The order parameters A and A are determined by the 
former equations, Eqs. (|6l7p . 

Thus the Gaussian variational approach is entirely con- 
sistent with the EOM method. However, it goes also be- 
yond it, since it enables us to obtain an estimate for the 
free energy. By Eqs. (|22[) and (j2"3")l . to calculate the best 
approximation Fq to the free energy, it is sufficient to 
consider local actions, for which we can express Sp, and 
thus Fq , in terms of a Hamiltonian 



H 



v 



5/A 



d> 



Ho -A 

2A+ 



2A 

(Ho - A* 



$ 



(24) 



Since the functional integrals are, by definition, normal 
ordered, the Hamiltonian H-p also needs to be normal 
ordered, as emphasized by the semi-colons in Eq. (|24| . 
indicating normal ordering with respect to the vacuum^ 
In this Hamiltonian language, Eq. (fTBf takes on the 
form 

F G (A, A) = -T log Zv + (H-Hv)-d, (25) 
with H the full Hamiltonian of the system, Eq. (fTJ) , and 

(26) 



Z v = Tre- /3ffD , 



= Tri 



,-PH- 



")/Z 



v ■ 



(27) 



Notice that Fq(A, A) also depends implicitly on the 
chemical potentials /i Q and the temperature T, and it 
must be minimized to find the mean field value of the 
variational parameters, A(/i Q ,T) and A.(fi a ,T). 

In this Hamiltonian approach, the evaluation of Eq. 
(1231) is straightforward (see Appendix [D]), and for the 
free energy density we obtain 



/ G =i|d£p(£)Tr(£-A) 



T 



d£p(£)Trlog(2cosh(/3B(£)/2)) 



(28) 



22 ((-A-a/3 - fJ. a S a p)n a p + A Q /3 (|n Q /3| 2 — n aa npp)) 



aP 



■ Y^ ( A "/9 d a/3 + A a/3 d "/3 " A Q/3 |d, 



ap\ 



aP 



Here j3 = l/T is the inverse temperature, the densities n 
and d arc determined by Eq. (fTO)) . and the matrix B(£) 
is defined in Eq. ([5]). 
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FIG. 5. (Color online) Comparison of the EOM and the 
Gaussian variational methods. Left: At low temperatures 
the EOM equations have multiple solutions, and become un- 
reliable close to first order phase boundaries. Right: The 
variational approach combined with a simulated annealing 
identifies correctly the physically relevant absolute minima 
of the free energy density, Eq. (|28|) . Parameters used were: 
A Q#3 PF = 0.112, 7 W = 0.5, T c /W = 0.011, and £ F /W = 
0.24 (half-filling). 



As stated before, at the local minima of the func- 
tional fa, the order parameters A and A fulfill the EOM 
self-consistency equations. In our numerical calculations, 
however, we have not enforced this constraint. Rather, 
we treated the order parameters as independent and free 
variables, and used a Monte Carlo method to find the ab- 
solute minimum of Eq. (|28p in the 15-dimcnsional space 
spanned by these order parameters. In the end, we ver- 
ified numerically that at the minima A and A indeed 
satisfy the EOM self-consistency equations. 

A comparison of the variational Monte Carlo ap- 
proach and the straightforward solution of the EOM self- 
consistency equations is presented in Fig. [S] At low tem- 
peratures, the EOM becomes unreliable in the vicinity 
of first order phase boundaries, and finds several possi- 
ble local minima. The variational Monte Carlo method 
(with simulated annealing), however, finds the absolute 
minimum of the free energy, fa, and is able to identify 
the physically relevant solution. 



C. Symmetries 

For an SU(3) symmetrical interaction, X a ^p = A, 
the structure of the phase diagram is largely determined 
by the underlying SU(3) symmetry. In particular, for 
fi a = /i the Hamiltonian is invariant under global SU(3) 
rotations, ^ a (x) H> J2p U a p^p(x), and a global U(l) 
gauge transformation, ^ a (x) i-> e ltp ^> a (x). 

The ferromagnetic order parameters n and A are Her- 
mitian. They are invariant under the U(l) gauge trans- 
formation, and transform under SU(3) rotations as 



n T h->Un T U t , A^UAU* 



(29) 



which, — after taking out the trivial trace, — is equiva- 
lent to the 8-dimensional adjoint representation of SU(3). 
The order parameters d and A arc, on the other 
hand, skew-symmetric, transform as d 4 e l2lp d and 
A M> e l2lp A under U(l) gauge transformations, and the 
global SU(3) group transforms them according to 



d^UdU T , A^UAU T , 



(30) 



which is equivalent to the conjugate representation of 
SU(3). This can be seen by introducing the 3 component 
vectors d Q = \ J2p 7 ea/3 T d/3 7 and A = ± £\ £ a/3 7 A/3 7 
by means of the completely antisymmetric Levi-Civita 
symbol e Q/ 3 7 . In this form Eq. (|30[) reads 



d4 U*d, AkU'A 



(31) 



In the special case, Xa^p = A and /i a = //, symmetry 
implies that the Ginzburg-Landau functional must be in- 
variant under the transformations (|29p and ([3"U]) , and the 
U(l) gauge transformation. The onset of superfluidity, 
however, spontaneously breaks the SU(3) <g> 17(1) sym- 
metry down to SU(2) (g) C/(l). This spontaneous sym- 
metry breaking is accompanied by the emergence of five 
Goldstonc modest 

The presence of the chemical potentials, jj, a p = S a p/j, a , 
obviously breaks the SU(3) symmetry. However, one has 
strong symmetry-dictated constraints on the Ginzburg- 
Landau functional even in this case, and the latter 
must be invariant with respect to the transformations in 
Eqs. ([29| and (|30|) . provided that fj, is also transformed 
accordingly, p, M> U/iU^ (see also Section Hv|) . In addi- 
tion, even in the presence of chemical potential differ- 
ences, SU(3) symmetry implies Ward identities^! relat- 
ing four-point expectation values and the ferromagnetic 
order parameter n as 

(fia ~ ftp) n a p=Y J ^ (A/J 7 - -W (* 7 *i*/3* 7 ) . 

7 

(32) 
From this identity (derived in Appendix [AJ it follows 
that n is diagonal for an SU(3) symmetric interaction. 
We remark that a similar approximate Ward identity can 
be derived within the Gaussian variational method (see 
Appendix [Bj , leading to the same conclusions. 

The off-diagonal elements of the chemical potential 
tensor, (1 describe tunneling between different hypcrfinc 
components, and they typically vanish in practical situa- 
tions. Under these restrictions, allowed SU(3) rotations 
generate essentially only permutations of the hypcrfinc 
labels, a, and the corresponding chemical potentials, \i a . 
On the (fj, x ,fi y ) plane, these permutations translate to 
C3 rotations and reflections, and give a two-dimensional 
representation of the S3 ~ C^ v group, implying a trian- 
gular symmetry of the phase diagram in this plane (see 
Fig El). 

In addition to the symmetries discussed above, for an 
SU(3) symmetrical Hamiltonian, the mean field equa- 
tions also have a certain particle-hole symmetry if the 



single particle density of states obeys g(£) = g(— £), and 
the chemical potentials are set to a value, fi —> /^haif , such 
that g is exactly half-filled. Under these conditions we 
can show (see Appendix [E| that the mean field solutions 
are symmetrical in the sense that for Sfj, a = fx a — /^haif 
and for Sfi a — > —S/j, a the superfluid and magnetic sym- 
metries are broken in the same channels and the order 
parameters are also equal apart from signs, global gauge 
transformations, and conjugation. In this special case, 
due to the additional permutational symmetry discussed 
above, the phase diagram exhibits a sixfold Cq v symme- 
try in the (5[i x , Sfiy) plane for traceless chemical potential 
shifts, 5fj,\ + <5^2 + <fy*3 = 0, (see Fig. [TJ) 

This particle-hole symmetry also emerges at the level 
of the Hamiltonian in certain cases. The half- filled at- 
tractive three component Hubbard model on a bipartite 
lattice 

h = -tJ212( a l« a i<* + h - c -) - -^ E£ n « - 2 )2 ' 

a (ij) i a 

e.g., has an exact particle-hole symmetry: it is invariant 
under the unitary transformation ai a -f-> sign(i) a ia , with 
sign(i) taking values ± for the two sublattices. Just as 
the mean field symmetry discussed in the previous para- 
graph, this exact symmetry relates the order parameters 
of the symmetry broken phases for ±Sfi a . We remark 
that, on a lattice, for stronger couplings, in addition to 
the SF/magnctic phases discussed here, other non-trivial 
phases may emerge (eg. charge density waves or trionic 
phases)] 16 ' 25 

Although the particle-hole symmetry discussed here 
holds only for a single and special chemical potential 
value, we found that for T c <C W higher order terms 
in the Ginzburg-Landau action are only sensitive to the 
immediate vicinity of the Fermi surface. As a result, 
particle-hole symmetry becomes an approximate symme- 
try with a good accuracy, whenever the slope of the single 
particle density of states vanishes, 7=0. For \x a = (x, 
A a ^/3 = A, and 7 = we thus recover a phase diagram of 
hexagonal symmetry within our numerical accuracy (see 
Fig. P. 



III. MEAN-FIELD PHASE DIAGRAM 

Let us now present the phase diagrams in the weak 
coupling limit, T c <C W, as obtained numerically, by the 
EOM and Monte Carlo methods presented in Section [Til 



A. Constant density of states (7 = 0) 

As we argued in the Introduction, except for the SU(3) 
symmetric point, a system of constant DOS always favors 
the formation of a SF phase in one of the pairing chan- 
nels (12), (23) and (31), having the smallest chemical 
potential difference. If the chemical potential difference 



T = T, 



7* = 0.907V 




-6-4-2 2 4 6 -C -i -2021 
ihjT c Hx/T c 



(23) 
(31) 
(12) 



a 



(i 



FIG. 6. (Color online) Mean-field phase diagrams at constant 
DOS, 7 = 0. Different SF phases are separated by first or- 
der lines. At T — 0.25T C SF-N transitions are of first order, 
whereas the y beco me of second order for T > f Sarma « 0.48T C 
(see Section IIIIC|I . Absolute values of components of the or- 
der parameter A are given in units of T c (see color code). 
Parameters at the SU(3) symmetric point: \ a jtpPF = 0.1, 
T c /W = 0.0076, and 7 W = &/W = 0. 



between the components forming the SF state exceeds 
a certain limit (known as the Clogston limits at zero 
temperature in case of two fcrmionic components), the 
system goes into the normal phase. This transition can 
either be of first or of second order, depending on the 
temperature^ 

Fig. [6] shows the numerically obtained phase diagram 
at different temperatures. All these cuts have the struc- 
ture presented in Fig. [TJ The hexagonal symmetry of the 
middle of the phase diagram is related to SU(3) symme- 



try: it is due to the invariance of the Hamiltonian un- 
der the permutations of the fermion species (a ■<->• /3 and 
fi a <H- fip) and the approximate particle- hole symmetry, 
as explained in Section lll CI The first order SF-SF transi- 
tions appear along lines where the chemical potential dif- 
ferences between two different pairs of fermions become 
equal. Along some special directions in the (/2 X , /j, y ) plane 
two out of three fermions have equal chemical potentials, 
and can form a SF state even far away from the central 
SU(3) symmetric point. This explains the ray-like struc- 
tures in Fig. [5] In all other directions the chemical po- 
tential differences continue to grow until the system goes 
into the normal phase at chemical potential differences of 
the order of the superfluid gap at the SU(3) symmetric 
point. For T > f Sarma w 0.48 T c this chemical potential 
driven SF-normal transition is of second order, however 
it becomes of first order below T' Sarma (see Section llll CI) . 



B. Linear density of states (7 7^ 0) 

In case of a non-constant DOS (7 ^ 0), particle-hole 
symmetry is broken at the Fermi surface, even at the 
SU(3) symmetric point. At a first glance, the phase 
diagram is only slightly different from the 7 = case, 
however, at a closer look qualitative differences can be 
discovered (see bottom and top parts of Fig. [7]). For 
7^0 the SF state not necessarily forms in the channel 
with the smallest chemical potential difference. The rea- 
son is that the gap is exponentially sensitive to the DOS. 
As a result, it may be favorable to form an SF state in 
channels, where the DOS is larger at the chemical po- 
tential, even at the expense of Zeeman energy (chemical 
potential) loss. This mechanism is driven by the deriva- 
tive of the DOS 7, and changes the phase diagram close 
to the SU(3) symmetric point. Here the phase diagram 
has only three-fold symmetry, corresponding to 'color' 
permutations, and superfluidity forms in channels of the 
largest density of states. At higher values of the chemical 
potential, however, the phase diagram remains essentially 
unaltered, and is almost identical to that of constant den- 
sity of states. 

These results are similar to the predictions of Ref. |24| . 
however, the phase structure differs somewhat, and the 
direction of the phase diagram of Ref. [24| seems to be 
flipped. We verified, that both the variational calculation 
and the equation of motion method yield consistently 
the phase diagram presented here, which we can also re- 
produce by the Ginzburg-Landau approach, presented in 
Section [TV] As we discuss there, the Ginzburg-Landau 
action of Ref. [2J cannot produce the six-fold symmetric 
structure of the overall phase diagram, and one needs to 
keep higher order terms to recover it. 

The previously discussed region of three-fold symmetry 
is, however, usually small compared to the overall scale 
of the phase diagram. For the parameters of the left 
figures in Fig. [7J e.g., T c /W = 0.011, and a relatively 
steep density of states with 7VF = 0.5, the three-fold 
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FIG. 7. (Color online) Phase diagrams at constant (left) and 
linear (right) DOS at T — 0.5T C . Non-zero 7 deforms the 
middle of the phase diagram (top right), whereas on the large 
scale (middle), the phase diagram with constant and linear 
DOS are almost identical. Parameters at the 5(7(3) symmet- 
ric point: (\ a ^3 P F, T c /W, jW,£ F /W) = (0.1,0.0076,0,0) in 
the left and (0.112, 0.011, 0.5, 0.24) in the right figures. 
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FIG. 8. (Color online) Phase diagrams at T — T c with 
constant (bottom left) and linear (jyW = 0.5) DOS (top 
left and bottom right). For linear DOS the SF-N critical 
temperature can exceed T c of the SU(3) symmetric point, 
whereas for 7 = 0, the SF phase disappears everywhere 
above T c . The largest values of the color scales correspond 
to \A a p\ = 0.07T C (top left), and \A a0 \ - T c (bottom left 
and right). Parameters at the SU(3) symmetric point are 
(Ac^pf, T c /W, jW,fr/W) = (0.1,0.0076,0,0) in the bot- 
tom left and (0.112, 0.011, 0.5, 0.24) in the top left and bottom 
right figures. 



symmetric region is present only for \fi x \, \(j, y \ < 0.1 T c , 
while the overall scale of the phase diagram is about ~ 
3 T c . The relative size of this central region increases for 
larger interaction strengths, and for T c /W = 0.105 and 
^W = 0.5 we find, e.g., that the central triangular region 
extends to |/i x |, \n y \ < 0.25 T c . The size of the central 
triangular region seems to scale roughly as ~ \FfT~c- 

In Fig. [8] we confirm the predictions of Ref. [24|, that 
breaking the SU(3) symmetry by the chemical poten- 
tial can indeed lead to the appearance of superfluidity. 
Again, this is simply related to the fact, that the su- 
perfluid transition temperature is exponentially sensitive 
to the DOS at the Fermi energy. At the SU(3) critical 
temperature T c , superfluidity appears only in small re- 
gions of the phase diagram, around the lines where two of 
the three fermion species have equal chemical potentials. 
These regions lie on that side of the SU(3) symmetric 
point, where the particles of the closest chemical poten- 
tials have higher DOS at the Fermi energy than the third 



one. We remark that the expansion of the free energy up 
to third order in the order parameters can not recover 
this structure precisely, and here the phase diagram is 
significantly different from the phase diagram of Ref. 1241 . 



C. Two component superfluidity 

It is instructive to compare our mean-field theory with 
results obtained for two component systems. As first no- 
ticed by Sarma^ for two component systems the Zecnian 
field-induced SF-N transition becomes of first order be- 
low the temperature T Sarma , and above the chemical po- 
tential difference ^ arma = (^f arma -//f arma )/V^. Sarma 
also determined the mean- field values of this critical point 
(Sarma point), and obtained 



nSarma 



0.58 T c 



Sarma 



1.5 T c 



(33) 



10 



with T c the critical temperature at fi x = 0. He also 
determined the critical chemical potential difference at 
zero temperature, known as the Clogston limits, 



^ciog = 2 A(T = 0) = 1.764 T c 



(34) 



with A denoting the SF order parameter. 

The three component system exhibits a two component 
behavior in regimes where the chemical potential of two 
species remains close, e.g. \ni— [Pz\ ~ T c , while that of the 
third component is very far from them (\ps\ ^> T c ). To 
investigate this limit, we fixed p y = 5T C , and varied p Xl 
along the solid line shown in the top left panel of Fig. [9] 
The corresponding SF phase diagram displays features 
similar to those predicted by Sarma. At T = tem- 
perature, the absolute value of the SF order parameter is 
independent of p x in the supcrfluid phase, and its magni- 
tude agrees with the BCS result, A(T = 0) = 0.882 T^*\ 

with Tc being the critical temperature at \x x = 04£ The 
critical value of p x (Clogston limit), however, shows sig- 
nificant deviations compared to Eq. (|34[) . For a coupling 
A = Xpp =0.1, e.g., we find both for a two and for a 
three component system 



/'• 



Clog 



"> #?° S |a=o.i = 2-19 A(T = 0) = 1.93T« . (35) 



For T c « f, the prefactor was found to be approxi- 
mately independent of the value of p y and particle-hole 
symmetry breaking parameter, 7. The difference be- 
tween Eq. ([35)) and Clogston's result is due to the in- 
clusion of magnetic degrees of freedom in the free energy 
density, Eq. (|2"5j). which accounts for interaction- related 
contributions to the Pauli susceptibility, \ ~ pp, ne- 
glected in Sarma's work^ These susceptibility contri- 
butions are proportional to \ a /3p F , and therefore result 
in a correction to the magnetic energy of relative size 
~ A12/OF, in rough agreement with the numerically ob- 
served shift of p x lo& - It is easy to understand this dif- 
ference on physical grounds: In the SF state (12), the 
densities nn and 7122 are exactly equal at T = 0, while 
in the normal state they shift according to the chemical 
potential difference. The interaction is, however, repul- 
sive in the magnetic channel. Consequently, the (mag- 
netized) normal state becomes less favorable, and p x los 
shifts upwards. 

Locating numerically the Sarma point we also find that 
it is shifted compared to Eq. (f3lt)) . 



T-iSarma , rp Sarma 

A=0.1 

Sarma . ~ Sarmal 

H'x ~* H'x Ia=0.1 



0.48 T" 



1.842 T, 



(*) 



(36) 

(37) 



again, approximately independently from the value of 7. 
These results and Eq. (|35p demonstrate that the posi- 
tions of the Sarma point and the Clogston point, Eq. (jM)) 
can significantly deviate from their standard BCS values 
due to interaction effects. Furthermore, their indepen- 
dence from the particular value of 7 shows that, at least 
for T c <C W, particle-hole symmetry breaking does not 




FIG. 9. (Color online) SF phase diagram (top right and bot- 
tom) at linear DOS (jW = 0.5), with p y = 5T C kept constant, 
as indicated by the solid line in the top left figure. The SF-N 
transition becomes from second order (solid line) to first or- 
der (dashed line) below the temperature f Sarma = 0.48 T C W , 
and chemical potential difference p^ arma = 1.842 T c , with 
Tc — 1.027 T c the critical temperature at p x = and 
p y — 5T C . Parameters at the SU(3) symmetric point were: 
X^ppp = 0.1, jW = 0.5, T c /W = 0.0076, £ F = 0. 



have a significant effect on the SF phases in the regime 
where the chemical potentials arc far from the SU(3) 
symmetric point. 

In the SF state, the SF species are bound together, 
and the condensate itself cannot be polarized. This has 
an experimentally important manifestation at the SF-N 
transition, where a sudden shift appears in the densities 
at the phase boundary, as presented in Fig. [TO] At zero 
temperature, the densities in the SF channel are equal, 
and their value does not depend on the chemical potential 
difference, whereas at the SF-N transition, a difference in 
the densities sets in. At temperatures below T Sarma 7 the 
SF-N transition is of first order, and the densities jump 
discontinuously on the phase boundary. In Fig. [TUJ this 
amounts toa~ 1% jump in the densities. In the strongly 
interacting regime, however, the jump is expected to take 
much higher values, similar to two component systems^ 

Let us close this section by investigating the effect of 
SF transition on the third, normal component. Indeed, 
in the presence particle-hole symmetry breaking, the SF 
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order parameter couples directly to the magnetization, 
and should shift the density of the third component. 
Fig. [10] shows this effect for a linear DOS in the weak 
coupling limit. We find that the shift in the density of 
the third component is only of the order of 0.01% for 
T c /W = 0.0076, however, for larger ratios, T c /W = 0.1 
(but the same 7) it reaches values of the order of 1%, in- 
dicating that this effect may be measurable in the strong 
coupling regime. 
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FIG. 10. (Color online) Interplay between superfluidity and 
magnetism in the SF channels, a = 1,2 (top), and for the 
third, normal component (bottom), for linear DOS (fW = 
0.5), with fx y — 5T C kept constant. The shift of the densities 
along the SF-N phase boundary is smooth for T > f Sarma 
(solid line), and discontinuous for T < T Sarma (dashed line). 
The density jump of the third component is much smaller 
than that of the SF components. [See also Fig. [9]] Param- 
eters at the 5(7(3) symmetric point: \ ai tpPF = 0.1, 7W = 
0.5, T c /W = 0.0076, & = 0. 



IV. GINZBURG-LANDAU ACTION 

In this section, we focus on the central region of the 
phase diagram, and construct a Ginzburg-Landau (GL) 
expansion of the free energy (|28[) around the 517(3)- 
symmetric point, fj, x = fi y = for T sa T c . Through- 
out this section, we assume a perfectly SU(3) symmet- 
rical interaction, X a ^p = A. While the form of the 
Ginzburg-Landau functional is dictated by symmetry, 
the coefficients of the various terms depend on the mi- 
croscopic parameters. We shall give approximate expres- 
sions for them, as obtained through a numerical analysis 
of Eq. (j2g|). 

In the weak coupling limit, the dimensionlcss free en- 
ergy density, 

Jg = !gI{pfTI) , 

can only depend on a few dimensionless physical pa- 
rameters: the dimensionless interaction A = pf\ the 
dimensionlcss slope of the DOS at the Fermi energy 
7 = jT c , the reduced temperature t = (T — T c )/T c , 
and the dimensionless chemical potential differences 6jl = 
(fi - p su ^)/T c , with ^ sc/ ( 3 ) denoting the chemical po- 
tential at the SU(3) symmetric point.— Most impor- 
tantly, however, Jq is a functional of the dimensionlcss 
order parameters, 

A = A/T c , 8 A = (A - A su{3) )/T c , (38) 

with A ' 3 ' denoting the renormalized chemical poten- 
tial at the SU(3) symmetric point. 

The expansion of the free energy contains only SU(3) 
invariant terms and can therefore be expanded as^ 

f G = A, Tr(AA+) + A 2 Tr((AA+) 2 ) (39) 

+ Bi Tr(o~A 2 ) + B 2 Tr(,5A) 2 + B 3 Tr(SJi SA) 
+ d Tr((5AAA+) + C* 2 Tr((5A)Tr(AA+) 
+ C 3 Tt(SJ1AA + ) + .... 

The 8 coefficients appearing in this expansion are all 
functions of A, i, and 7. We determined them by fit- 
ting the free energy Eq. (|28)) numerically, and found that 
the expressions in Table UJ give a good estimate for these 
parameters 42 At the minima of the free energy functional 
above we have SA oc SJ1 and A oc \ft. Therefore, the ex- 
pansion above contains all terms up to 0{t 2 , SJlt, SJ1 2 ) . 

The superfluid phase transition is driven by the term, 
Ai(t, A), which changes sign at the SU(3) point. All 
other coefficients are approximately constant close to the 
phase transition. The terms ~ Bi describe the ferromag- 
netic order parameter, and its response to the external 
" magnetic field" , ju. The most interesting terms are the 
ones proportional to the coefficients C^: these describe 
the coupling between the SF order parameter and the 
magnetization (or chemical potential differences), and 
they are responsible for the three-fold symmetric struc- 
ture in the central region of the phase diagram (see 
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parameter 


approximate expression 


A 1 


2.00 t + ... 


A 2 


0.40-1.20*+... 


Si 


0.5000 + 1.000 A +... 


B 2 


-1.000 A + ... 


B a 


-1.000 + ... 


Ci 


1.257 + .. . 


c 2 


-1.227+.. . 


c 3 


-0.62 7/A + ... 



TABLE I. Approximate expressions of the Ginzburg-Landau 
coefficients in Eq. (J39J). The dimensionless parameters are 
A = p F \, 7 = 7 T C , and t = (T - T c )/T c . 



Fig. [7]). The terms C\ and C2 couple the superfluid and 
magnetic order parameters, and produce the density shift 
of the normal component at the onset of superfluidity. 
Notice that all these terms are found to be proportional 
to the dimensionless particle-hole symmetry breaking pa- 
rameter, 7. 

While the third order expansion, (|39p accounts for 
the central regions on the right panels of Fig. it 
does not recover the sixfold symmetric structure that 
dominates the phase diagram at larger chemical poten- 
tial differences. This is obvious, since the terms C\,C<i 
and C3 are odd under the particle-hole transformation, 
SJI <-> —SJI, SA -f-> — <5A*, and are proportional to 7, while 
the hexagonal structure is even under particle-hole trans- 
formation, and already appears for 7 = 0. The "hexag- 
onal" structure must therefore be controlled by higher 
order terms, containing even degree polynomials of SJI 
and SA, coupled to the SF order parameter. Unfortu- 
nately, the number of such terms is huge, and is next to 
impossible to determine all of them and their correspond- 
ing GL coefficients accurately. However, observing that 
the ferromagnetic response is always small, we can just 
focus on the SF order parameter. At a formal level, this 
can be done by minimizing the free energy functional fc 
in 5 A for any fixed SJI and A, and thus defining 

f G (SJl,SA) = f G (SJi,SA,SA min (5Ji,SA)) . 

The form of this GL functional is also dictated by sym- 
metry, and it can also be expanded in SJI and SA. Up to 
0{t 2 ,5Jl 2 t) itrcadsi 2 . 

f G = 01 Tr(AA+) + a 2 Tr((AA+) 2 ) 

+ 6Tr(^AA + ) + ciTr((5/I 2 AA + ) (40) 

+ c 2 Tr(SJlASTiA + ) + .... 

The approximate values of the numerically obtained co- 
efficients are enumerated in Table ILT1 

Minimization of Eq. (|4T))) yields the correct structure of 
the phase diagram in the vicinity of the SU(3) symmetric 



parameter 


approximate expression 


ai 


2.0 t + ... 


a2 


0.40- 1.2t+ ... 


b 


(3.2i-0.083/A 2 ) 7 + ... 


C\ 


0.125 -0.29 A- 0.13 1 + ... 


C2 


-0.115 + 0.27 A + 0.12 £ + ... 



TABLE II. (Color online) Approximate expressions of the 
Ginzburg-Landau coefficients in Eq. (|40|) . The dimensionless 
parameters are A = pfA, 7 = *yT c , and t = (T — T c )/T c . 



point, and accounts for the competition between the odd 
(&,...) and even (ci,C2, ...) order couplings. We also 
checked that it determines correctly the absolute value 
of the SF order parameter in the weak coupling regime 
T c /W < 0.1 at temperatures 0.9T C < T < T c . However, 
the locations of the triple points at the interface of the 
threefold and approximately sixfold symmetric structures 
in Fig.UJare reproduced only with an error of about 50%. 
Although this error is very large, it is also natural, since 
on the scale of this structure, the chemical potential dif- 
ference is of the order of SJI rs 0.2. Therefore SJI cannot 
be considered as a small parameter here, and higher order 
terms in the expansion (|40[) shift the phase boundaries 
significantly. 



V. BEYOND MEAN-FIELD 

In the discussion presented so far we restricted our- 
selves to a mean-field approach, and neglected fluctu- 
ations. Fluctuations, however, not only reduce some- 
what the transition temperatures and fields, but they 
also change the universality class and thus the critical 
exponents of the transition. In ordinary superfluids, such 
fluctuation effects are typically hard to observe, however, 
in cold atomic systems one can reach the strong coupling 
regime, and therefore a non-trivial critical behavior may 
be observable^ 

First, let us discuss the central SU(3) symmetrical 
point of the phase diagram, fj. x = \x y = 0. At this point 
only the first two terms of the GL action (|4"D)) survive for 
an SU(3) symmetrical interaction. These terms as well 
as the gradient term, Tr{d r A ■ d r A + } have an increased 
0(6) symmetry with respect to SU(3)^- with the real 
and imaginary parts of the independent components of 
A forming a six component real vector. Since higher 
order terms are irrelevant in the renormalization group 
(RG) sense, the fi x = [i y — transition is described by 
the 0(6) critical theory. Thus the correlation length di- 
verges as £ ~ \T — r c | _ly °( 6 ), while the order parameter 
scales as (A) ~ \T - T c f°m. For d = 3 dimensions, 
the critical exponents are known from e expansions, 45 
\/n expansions^ as well as from high-temperature 
expansions^ and Monte- Carlo simulations js& giving sim- 
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ilar results, 



,,3D 

y o(6) 



0.80, 



A 3D 



0(6) 



0.41 



In two dimensions, on the other hand, fluctuations sup- 
press the phase transition at the SU(3)-symmetrical 
point, T 2D — > 0j££ which thus becomes a quantum criti- 
cal point. 

For generic values of fJ. x ,fJ-y ^ 0, only one superfluid 
channel dominates the phase transition, which is there- 
fore described by the XY model. In d = 3 dimension the 
corresponding critical exponents are given by^ 



3D 
'XY 



0.67, (3 



3D 

XY 



0.35, 



(41) 



while in d = 2 dimension the transition is of Kostcrlitz- 
Thoulcss type. 50 

Interesting critical behavior emerges in the vicinity 
of the bicritical lines of Fig. [2] Along these lines, two 
components of the matrix A, e.g. A13 and A23 com- 
pete with each-other to form the superfluid. These can 
be grouped into a real four component vector, p = 
(ReA23, ImA23, ReAi3, ImAi3). Fermion number con- 
servation implies that the effective action must be invari- 
ant under global phase transformations, A,-j — ► e'^Ay, 
which translates to an 0(2) x 0(2) symmetry in terms of 
the field ip. Up to fourth order, the most general effective 
Hamiltonian can be written as^ 

H L w= Jd d x [i(V^) 2 + t+p 2 +t^pUp (42) 

+ u{p 2 ) 2 + v(pUp) 2 + wip 2 (p>np) + ...], 

where the terms breaking the 0(4) symmetry were writ- 
ten in terms of the matrix 



n 



(1 



\ 



-1 



(43) 



-v 



In the absence of the terms t— and w, this action has 
an additional Z2 symmetry, A 13 <-> A23, leading to a 
0(2, 2) = (0(2) x 0(2)) x Z 2 symmetry of the free energy 
functional. In the presence of particle-hole symmetry, 
one can show that at the boundary of the two superfluid 
phases the Z2 violating terms vanish: t— = w = 0. In 
general, however, the simultaneous vanishing of i_ and 
to is not guaranteed. Nevertheless, already leading order 
e expansion indicates^ that the coupling w is irrelevant 
at the phase transition, t± — > 0. Thus the Z2 symmetry 
is apparently restored at the transition, and the critical 
state must be described by the 0(2, 2) symmetrical func- 
tional with t-,w —>■ 0. 

The 0(2,2) functional (02) with t-,v) -> thus de- 
scribes the phase transition at all bicritical endpoints 
where two superfluid phases meet (white circles in 
Fig. [TTj) . Notice that the structure of the phase diagram 
changes close to T c , and the six 0(2, 2) points, - char- 
acteristic at lower temperatures, - pairwise merge into 



T C >T> T tn T tTl > T 



(23) 




31) 




f^/T e 



FIG. 11. (Color online) Schematic picture of the position of 
the 0(2, 2) points (empty circles) in case of linear DOS. At the 
temperature T trl below which the triple points appear, from 
each 0(2, 2) bicritical line (left) two new bicritical lines of 
the same universality class branch out (right). The branching 
points are multicritical. SF-SF transitions are of first (dashed 
lines), whereas SF-N transitions are of second order (solid 
lines). 



three 0(2, 2) points above a tricritical temperature, T tn , 
as also shown in Fig. ITT]) . 

The second order terms t+ and £_ trigger the SF-N 
and SF-SF transitions, respectively, and scale as 



oc (5/X||, 
oc<5^_i_, 



(44) 
(45) 



for small chemical potential shifts parallel (5/a\\) and per- 
pendicular (Sfi±) to the SF-SF phase boundary. 

The model (|4"2")l has been studied extensively ; 51 ' 52 typ- 
ically in the framework of the more general n ■ m com- 
ponent models^ Despite the extensive effort, the stabil- 
ity of its various fixed points is still debated. System- 
atic e expansion yields three non-trivial fixed points with 
t*_ = w* = 0, which could potentially describe the criti- 
cal state: (a) an 0(4) Heisenberg fixed point with u* > 
and v* = w* = (b) a decoupled fixed point (DFP) 
(u* = v*,w* =0), where the two superfluid components 
are described by two independent XY theories, and (c) 
a mixed (or biconical) fixed point (MFP) with u* ^ v* 
and w* = 0. 

For small values of e = 4 — d, e expansion yields the 
picture shown in Fig. 1121 predicting that the mixed fixed 
point (MFP) describes the phase transition along the 
0(2, 2) critical line. However, already in second order 
in e^ the fixed point structure changes completely as 
one approaches the physical value, e = 1, and even the 
results of six loop e expansion remain completely incon- 
clusive regarding the stability of the fixed points^ 5 - Non- 
perturbative arguments, on the other hand, seem to sup- 
port that the rather boring decoupled fixed point (DFP) 
describes the critical stated 5 - - — 

The universality class of the fixed point has consid- 
erable impact on the phase diagram. The ratio of the 
critical exponents y± associated with the terms t± de- 
termine e.g. the shape of the SF-N phase boundary in 
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the vicinity of the bicritical point. Standard cross-over 
scaling arguments^ lead to the conclusion, e.g., that the 
specific heat diverges in the vicinity of the SF-N transi- 
tion line as 



Vt + < Vu 



c (t+)«|t + -C(t_)|- 



(46) 



where axy denotes the specific heat exponent of the XY 
model, and the phase boundary is determined by the 
function C(t~) 



C(t-)<x \t^\ v+/y - 



(47) 



Since the critical exponents y± arc different for the two 
possible stable fixed points even to first order in e, 



yf FP = 2-e/2- 



MFP 



V- 



= 2 - e/6 • 



(48) 



the shape of the phase boundary will be different in the 
two cases. Notice that since the DFP describes two in- 
dependent XY models, its exponents y± FP will be equal 
to all orders in e, implying that the SF-N boundaries 
start linearly at the bicritical point. For the MFP, on 
the other hand, y+ FP < y^ FF , and the SF-N boundary 
has a universal exponent in the vicinity of the 0(2, 2) 
point, as shown in Fig. [T31 This difference in the shape 
of the phase boundary provides a clear fingerprint of the 
universality class of the transition. 

The critical exponent /3 of the order parameter, (ip) 
along the SF-SF phase boundary, is determined by the 
exponent yh of the "magnetic field" at the critical fixed 
point, 



P 



Vh 



v+ 



Since the magnetic field exponents get their first non- 
trivial contribution in 0(e 2 ) order, to leading order in e 




FIG. 12. (Color online) Schematic picture of the 0(e 2 ) RG 
flows in the (u, v) plane, with t- = w = 0, for e < 5/7. For 
e — > 1 the fixed point structure changes, and the e expansion 
is inconclusive. 



y*+ = yt- 



T, (23) 




fix/Tc 



FIG. 13. (Color online) Schematic phase diagram of the vicin- 
ity of the 0(2,2) bicritical point (see Fig. [TJ in case of the 
mixed (left) and the decoupled (right) fixed point. In the for- 
mer case fluctuations modify the SF-N phase boundary into 
curves with universal scaling. 



we have 



Vh 



,MFP 



»"5 



(49) 



However, since y+ FP ^ 2/+ FP , the exponents /3 DFP and 
/3 MFP turn out to be different already to first order in e, 



P 



DFP 



1 3 

~ 2 ~ 20 £ 



/3 



MFP 



(50) 



VI. EXPERIMENTAL RELEVANCE 



Currently maybe 6 Li ultracold gases provide the most 
promising perspective for the realization of three com- 
ponent superfluidity. For high magnetic fields, the s- 
wave scattering lengths between the three lowest hyper- 
fine states approach the spin-triplet scattering length, 
ai2 ~ ^23 ~ 031 ~ — 2140 ao, with ao the Bohr radius^i 
At fields of ~ 2000 G, for example, the scattering lengths 
all deviate less than 2% from their average valued. It has 
been proposed theoretically that this deviation can fur- 
ther be decreased using radio frequency and microwave 
fields^ down to 0.1%, and thereby a strongly attrac- 
tive system can be realized with almost perfect SU(3) 
symmetry in this high field regime. 

Although three-body loss is a major obstacle in three 
component experiments, recent experiments showed that 
decay rates tend to decrease at high fields in 6 Li systems, 
and indeed, Fermi degeneracy has successfully been re- 
alized in this three component systemJi A 6 Li experi- 
ment on a system of Fermi energy Tp = 1/iK and with- 
out optical lattice would correspond to the parameters 
Ac^/3/Of w 0.11, 7W w 0.18 and T c /W « 0.01,2a This 
system would thus be in the regime of weak interactions, 
studied here. However, such a small critical tempera- 
ture is currently unreachable. Application of an opti- 
cal lattice can, however, easily bring the system into the 
regime of strong interactions, where SU(3) superfluidity 
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may be accessible. Though our calculations do not ap- 
ply for strong interactions, we believe that, similar to the 
SU(2) case,^£&22 the major features of our phase dia- 
gram are robust, and should carry over to the strongly 
interacting case. 



T = 0.97T, 



T = 0.907i 12) 




FIG. 14. (Color online) Phase diagram with only approx- 
imately equal interaction strengths. The rate of the criti- 
cal temperatures in the respective channels are T c /T c = 
T c (31) /Tc (12) = 0.95. The SF phase in the channel of strongest 
interaction repels the other two phases from the central re- 
gion of the phase diagram. Parameters at the fi x = fi y = 
point: (A 12 ,A 2 3,A 3 i)pf = (0.1057, 0.1046, 0.1046), T C (12) VF = 
0.01, 7W = 0, fr = 0. 



So far we assumed a perfectly SU(3) symmetrical in- 
teraction in our calculations. The phase diagram is, how- 
ever, somewhat modified if the the scattering lengths are 
only approximately equaL 61 ' 62 In Fig. [T?] we present a 
phase diagram for the case where we have set the ratio 
of critical temperatures in the different channels to be 



T c (23) /T c (12) 



(31) /T ,(12) 



Tr>/T ( 



0.95. For T^ ij) /W ~ 0.01 



this would correspond to a ~ 1% asymmetry of the scat- 
ring 

the SF phase is formed only in the (12) channel. 



(10) (231 

tcring lengths. At temperatures T c > T : > T c 
T ( 

The star-like shape of the phase diagram is preserved at 
lower temperatures, however, the interaction asymme- 
try destroys the sixfold symmetry of the central region 
of the phase diagram, including the 0(6) critical point: 
the phase (12) dominates this central region and expels 




FIG. 15. (Color online) Possible trap configurations for total 
atom numbers JV3 > N? > Ni. 



the other two SF phases. Thus the shape of this region 
depends rather sensitively on the interaction asymme- 
try, and fine tuning of the scattering lengths (by using 
RF and MW fields^ e.g.) may be needed to realize an 
SU(3) symmetric superfluid. 



VII. CONCLUSIONS 

In this paper, we studied the phase diagram and the 
interplay of fcrmionic and superfluid order parameters in 
a three component fcrmionic mixture. We mostly focused 
on the case of SU(3) symmetrical interactions, and stud- 
ied the weak coupling regime, where the critical tempera- 
ture is much smaller than the Fermi energy of the atoms, 
T c < Ep. We combined two complementary mean field 
methods (Gaussian variational method, and equation of 
motion techniques) to study how a chemical potential im- 
balance polarizes the atomic cloud and modifies/destroys 
superfluid order. Though the phase diagram of the three 
component system is naturally much richer than that of 
the two component mixture* 2 ^ there are some similari- 
ties: large chemical potential imbalances {\fJ>i — Hj\ 3> T c 
for all i y^ j), for example, destroy superfluid (SF) order, 
similar to two component mixtures. The corresponding 
SF-normal transition is of second order at higher temper- 
atures, while it becomes of first order below the Sarma 
temperature. 

The superfluid phase is, on the other hand, much richer 
than in the two component case. SF order can form in 
channels (12), (23), and (31), and the chemical potential 
driven transitions between these superfluid phases are of 
first order. In a real experiment, where fermion num- 
bers are approximately conserved for each component, 
such first order transitions would appear as segregation 
of different SF phases, and domain formation^ Experi- 
mentally, these domains would probably appear as a shell 
structure, sketched in Fig. \H>\ For N 3 > N 2 > N±, e.g., 
one expects that in the center of the trap a (32) super- 
fluid forms, however, approaching the external region of 
the trap T c decreases, and the (31) superfluid state be- 
comes more stable. 

As a rule of thumb. SF order tends to form in the 
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channel of the smallest chemical potential difference. 
This simple rule determines the overall structure of the 
phase diagram (see Fig [7]). However, unlike the two 
component case, for three component mixtures a non- 
trivial coupling between magnetic and SF order is also 
allowed . 16 ! 24 This interesting coupling — the strength of 
which is regulated by particle-hole symmetry breaking, 
7 ~ g'(Ep)/ ' g(Ep) ~ 1/Ep — leads to a peculiar trian- 
gular structure in the central region of the phase diagram, 
fii ~ /i, in agreement with the predictions of Ref. |24| 
(though with opposite orientation, see Fig [7]). However, 
the relative size of this central region is apparently pro- 
portional to ~ \/jT c ; therefore, for weak and interme- 
diate couplings, the triangular structure appears only in 
the close vicinity of the SU(3) symmetrical point, (U, = [i. 
For very strong attractive interactions, T c w Ep ~ W, 
on the other hand, the central (triangular) region must 
get more extended, and may become observable. 

We also constructed the Ginzburg-Landau functionals 
describing the three component mixture, and determined 
the temperature and asymmetry (7) dependence of the 
various coefficients. We have shown that, to capture the 
termination of the central triangular region, one needs 
to go beyond the expansion of Ref. [24|, and higher order 
terms need be incorporated in the functionals. 

As discussed in Sec. [V] fluctuations modify somewhat 
the mean-field picture. The temperature-driven phase 
transition for generic (unequal) chemical potential val- 
ues is typically described by the XY model and its critical 
exponents. However, for certain special chemical poten- 
tials, the competition between various superfluid orders 
may lead to interesting critical behavior. For pn = jj, and 
an SU(3) symmetrical interaction, e.g., the normal-SF 
transition belongs to the 0(6) universality class, and is 
characterized by the corresponding exponents. Along the 
critical lines separating the three phases, (12), (23), and 
(31), on the other hand, an interesting 0(2,2) critical 
behavior may emerge (see our discussion in Sec. IVp . The 
shape of the phase diagram in the vicinity of these special 
lines is then determined by the corresponding universal 
cross-over exponents. We emphasize that - while it is 
very difficult to observe it in the weak coupling regime 
- a non-trivial critical behavior could be observable in 
the strong coupling regime, often reached in cold atom 
experiments. 

Finally, we studied the fragility of the SU(3) physics, 
i.e., the sensitivity of these results and the phase diagram 
to the symmetry of interaction. We have shown that 
already a small difference in the scattering lengths can 
substantially distort the SU(3) phase diagram, and the 
SF phase of the channel with the strongest interaction 
may suppress and mask the SU(3) symmetrical (0(6)) 
critical regime. These results agree with those obtained 
in Ref. |62|. Here, however, in contrast to Ref. |62|, we 
focused on the consequences of SU(3) symmetry (rather 
than on the consequences of its violation), and the ef- 
fects of the coupling between ferromagnetic and super- 
fluid order parameters, neglected in Ref. |62j. In addi- 



tion, we also discussed the role of fluctuations and the 
structure of the emerging critical states and multicriti- 
cal lines. Our results as well as those of Ref. [62| indicate 
that in experimental realizations, to observe the SU(3) 
physics, one should use systems with almost perfectly 
symmetrical interactions, similar to Yb^., or one should 
use some tricks to make all scattering lengths equal as 
much as possible^ Moreover, one should possibly stay 
in the strong coupling regime, T c ~ W, where the impact 
of a small asymmetry in the interaction is not exponen- 
tially large. 
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Appendix A: Exact Ward identities 

In this Appendix, by making use of the global SU(3) 
invariance of the functional measure, we derive exact 
Ward identities that give constraints on the possible val- 
ues of the order parameters and densities, Eqs. (|6]|9]). 

Consider the partition function Z, defined in Eq. (|12l) . 
For the current calculation we rewrite the action 
Eqs. (|13ll4j) in the form 



£(>(£) =^ / dxifi a ((d T +Ho)8 a i3 -fr a p)i>p, (Al) 

a/3 J 

S'int(r) = - y^ T aPl& I dxil> a ippi/}^s • (A2) 

a/3"{5 

by introducing £i a p = fi a d a p and T a p 7 s = 
|A Q/ g (SasSpj — SajSps)- An SU(3) transformation of the 
fields ipa(x) — > J2b U a /3ipp(x) translates to the transfor- 
mation of jl and r in the functional integral. Expressing 
U = cxp(i Yl =1 T] a T a ) with the Gell-Mann matrices T a , 
we find 






I7»=0 



= «E 7 (a Q7 ^-t- 7 a 7/3 ),(A3) 



17" =0 



2i (A Q/3 - Xy S ) S aS T| . (A4) 



The invariance of the functional integral with respect 



to global SU(3) transformations, J^r 
the Ward identity 



= 0, leads to 



Va=0 



(Ma - M/8) ^7— - 2 2_^ (A / g 7 - A Q7 ) 



dpLafl 



OT 



70/^7 



(A5) 



for any a and /3, from which Eq. (|32[) follows. 
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Appendix B: Ward identities in the Gaussian 
approximation 

Here, we derive approximate Ward identities, similar 
to those in Appendix \K\ that hold in the Gaussian ap- 
proximation. As explained in Appendix [C] we can as- 
sume that the inverse propagator in the definition of the 
partition function Z?>, Eq. (fT5]) , is local, 



DipDi^e 



iJ-dl^lJD-^lMl) 



(Bl) 



where D ~ 1 is defined in Eq. 

An SU(3) transformation of the fields ip a {x) — > 
^aUapippix) translates to the transformation of order 
parameters, 



A^UAU+, 
A h^ UAU T , 



(B2) 
(B3) 



see Eqs. (|29l30p . Using the invariance of the partition 
function with respect to these global SU(3) transforma- 
tions, we get the following constraints on the densities, 



Tr 




= , (B4) 



with T a A = [A, T a ] and T^ = 2 (T Q A + AT a *).^i Here 
the matrices T a , a = 1, .... 8, arc the Gcll-Mann matri- 



ces. 



In case of SU(3) symmetric interactions, at the solu- 
tions of the EOM equations, Eqs. (|6|7I10|) . this equation 
simplifies to the same form as the exact Ward identity, 
Eq. ([SID, 



(p a - fifj)n af j = 



(B5) 



Therefore, when neither two of the chemical potentials 
are equal, the matrix of densities n and that of renor- 
malized chemical potentials A arc both diagonal (see 
Eq. ©). 



Appendix C: Saddle point equation in the Gaussian 
approximation 

In this Appendix, starting from the saddle point 
equation, Eq. (|2"Tj) . we derive the saddle point form 
of the propagator T> in the Gaussian approximation, 
Eqs. (J22I23J) . Wc will use the notations of Section iHBl 

First, wc fix the arbitrariness in the form of T)~ l in 
the definition of S v , Eq. (TT7|). We split P" 1 into 3x3 
matrices 



T>-\1,2) 



T A {xi,x 2 ) r B (xi,.T 2 ) 
r c (xi,X2) T D (x 1 ,x 2 ) 



(CI) 



It is easy to see, that because of the anticommutation of 
the fields tp a and ip a , modifications of T>~ 1 that leave 



r^(x 2 ,xi), T B {xi 1 x 2 ) - r^(x 2 ,xi) and 
r^(x2,Xi) invariant, will not change Sx>- 



T A {x 1 ,x 2 ) 

r c (xi,x 2 ) - 3. c 

Therefore wc may assume that T>~ 1 

symmetry 



has the symplcctic 



(V- 1 ) (x 2 ,Xl). 



(C2) 

The saddle point equation, Eq. (f2"Tj) , gives very strong 
constraints on the form of P. In particular, it is equiva- 
lent to the EOM self-consistency equation of Scction lll Al 
To see this, we use the definition Eq. (fl7il) to rewrite 
Eq. ((2TJ) in the form 



1 SZ V = S(S - S v ) q 
Zv 8V{\,2) - £D(1,2) 



(C3) 



The calculation of the left hand side of this equation is 
straightforward. Using only the definition of Zx> (see 
Eq. OH)), and Eq. flUJ}, we get 



SZ 



v 



Zt> <JP(1,2) 



-v-\%\) 



(C4) 



To evaluate the right hand side of Eq. (|C3[) . omitting a 
constant term, we can write 

(S - S v ) v = -^ J dld2V- 1 (l,2)V(2,l) + (S int ) v . 

(C5) 
Then, it is easy to see that, the saddle point equation is 
equivalent to 



D _1 (l,2) = £> ~ 1 (l,2)-2 



S(S; 



int I'D 



6V(2,1) 



(C6) 



Expanding (Si nt )v using Wick's theorem gives a product 
of equal time propagators, whose variation according to 
the propagator matrix T> can be straightforwardly calcu- 
lated. We get the desired formulas, Eqs. (|22l23p . with 
the order parameters A and A satisfying the EOM self- 
consistency equations, Eqs. (|6I7[) . and (|T0| . This means, 
that the EOM method is consistent with the Gaussian 
variational approach. 



Appendix D: Calculation of the Gaussian 
approximation to the free energy 

In the following we calculate the Gaussian approxima- 
tion of the free energy, Eq. (|2"5]) . We first introduce the 
Fourier components V'a(r) = -4? X)k elkl " a <*k' obeying 

the anti-commutation relations {o^, at/ /?} = ^a/3<^kk', 
where fi denotes the volume. With these, the Hamilto- 
nian, Eq. (|24|) . takes on the form 



Hx> = 




a£,a_ k B(£ k ) 



ak 



Tr (£ k 
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with B(£) defined in Eq. ([5]), and the last term originating 
from normal ordering. 

From the above form, the calculation of Z-d = 
fp re -/3i? D j g straightforward, though some care is needed 
to avoid double counting in momentum space. Note that, 
because of the symplectic symmetry, Eq. (fTTj) , and Her- 
miticity of the matrix B(£), its eigenvalues are real and 
come in pairs. To each eigenvalue ??(£) there is another 
eigenvalue — t/(£). Using this property, logZ-p simplifies 
to 



log Z v =\ Y, Tr log (2 cosh ( §B (&) 



(D2) 



E Tr (^-A). 



The calculation of (H — Hx>)v is also straightforward 
using Wick's theorem. One finds 



o 



aP 



+ 2J (-A-o/3 - PaSap) n a/ 3 



(D3) 



a/9 



■ 2J Agpdgp + A* a/3 d a t 



Thus, using Eqs. (|D2|D3|) . we get the result Eq. §28$) for 
the Gaussian approximation of the free energy density . 

In order to evaluate Eq. (|D3|) . the densities and anoma- 
lous densities, n and d, also have to be determined. 
These can be easily calculated from the variations of (|D2j) 
with respect to A and A, leading to the same equation, 
Eq. (fTO)) . as the EOM self-consistency equations. 



Appendix E: Particle-hole transformation 

Particle-hole symmetry introduces a Z2 symmetry of 
the mean-field phase diagram, when the band is half- 
filled, the DOS is particle-hole symmetric (p(£) = p(— £)), 
and the interaction has SU(3) symmetry (\ a ^p = A). 
This symmetry together with the permutation symmetry 
of the fermion species makes the phase diagram six-fold 
symmetric, see Fig. Q] 

In this Appendix we calculate the effect of the particle- 
hole transformation 



*„(&) <— ► & a {x) 



(El) 



on the order parameters A and A. This transformation 
leaves the interaction invariant, whereas it modifies the 
bare chemical potentials and the single particle energies 



as 



%i) 



—M01 



(E2) 

(E3) 



p a r fl a 4A ?iniax] 

where n max = J_ w d£ p(£) is the density of the com- 
pletely filled band. The bare chemical potentials remain 
unchanged on the mean-field level at 

Mhaif = -2An max , (E4) 

which is precisely the condition for the band being half- 
filled (sec Eq. ©). 

In order to investigate the inversion symmetry of the 
phase diagram, consider two Hamiltonians with opposite 
differences in bare chemical potentials from half-filling, 

H (1) = H(Ho, Mhalf + S/ia, A, *!, tf a ), (E5) 

H^ = H{Ho, Mhalf - Sn a , A, *L * a ), (E6) 

as defined in Eq. ([]}. A particle-hole transformation of 
H^> leads to the equation 



H& = H(-H , 



Mhalf 



6p a ,X,^ a ,^ a )=H^, (E7) 



where ^ a = $!' a . Accordingly, the densities in the origi- 
nal and the particle-hole transformed system can be con- 
nected as 



M = /St 



(2)* 



$M*i(s)*^)>(3) = -<r + 



,(3) 



(2)* 



d^:> = (* a {x)*p(x)) {3) = -d>:> 



(E8) 
(E9) 



Then, it is also straightforward to show from the defi- 
nitions Eqs. (|6l7p . that the relation between the order 
parameters are 



,(2) 



-A (3) * A (2) = -A (3) * 



(E10) 



Looking at their definitions, we see that the only dif- 
ference between H^ 2 ' and H^ is in the sign of 'Ho- How- 
ever, if the DOS is electron-hole symmetric, 



p(0 = p(-0, 



(Ell) 



then all of the EOM self-consistency equations 
Eqs. (|6I7I10|) , and the mean-field free energy Eqs. (|10I28|) 
are identical in the two systems. Therefore, the set of the 
possible mean-field configurations have to be the same 
(A« = A< 3 ), AW = A ( 3 )). Putting this, and Eq. (|E10|> 
together, we obtain the desired equations 

A(/i h alf + Spa) = -A*(jU h alf - 5p a ), (E12) 

A(/i h alf + Spa) = -A*(/i hal f - 6 p a ), (E13) 

connecting order parameters at opposite Sp a values, with 
the other parameters of the system unchanged. 

We remark, that in the special case when dpi + 8pi + 
dps = 0, the particle- hole symmetry connects the points 
of the same (p x , p y ) plane, and the mean- field phase di- 
agram has an inversion symmetry. Away from this plane 
the inversion symmetry is only approximate, due to loga- 
rithmic corrections to the values of the order parameters, 
coming from the asymmetric cut-off. 
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